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Tj- \ In order to assess the accuracy of commonly used approximate exchange-correlation 
^ ' density functionals, we present a comparison of accurate exchange and correlation poten- 
O ' tials, exchange energy densities and energy components with the corresponding approx- 
^ . imate quantities. Four systems are used as illustrative examples: the model system of 
two electrons in a harmonic potential and the He, Be and Ne atoms. A new ingredient 
\ in the paper is the separation of the exchange-correlation potential into exchange and 
correlation according to the density functional theory definition. 



1. INTRODUCTION 

Within density functional theory (DFT), the ground state energy of an interacting 
system of electrons in an external potential can be written as a functional of the ground 
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' state electronic density |l|| . When compared to conventional quantum chemistry methods, 
^ this approach is particularly appealing since it does not rely on the knowledge of the 
complete A'^electron wave function but only of the electronic density. Unfortunately, 
although the theory is in principle exact, the energy functional contains an unknown 
quantity, called the exchange-correlation energy, E^c [p] , that must be approximated in 
any practical implementation of the method. Several approximate exchange-correlation 
functionals have been proposed in the literature, the most commonly used ones being the 
local density approximation (LDA) and the generalized gradient approximation (GGA). 

The local density approximation is the simplest and most widely used exchange-cor- 
relation functional: 

^X°^[P]= /p(r)4°^(p(r))dr, (1) 

where ^^^^ip) is the exchange-correlation energy per particle of a homogeneous electron 
gas of density p. E^^^ [p] is by definition the exact functional for a homogeneous electron 
gas and has been shown to give also a qualitatively good description of the ground state 
properties of a variety of highly inhomogeneous systems [H]. However, LDA does not 
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always provide sufficiently accurate results. For example, it almost always overestimates 
the binding energy and underestimates the bond-length of weakly bound molecules and 
solids 0. LDA fails also to predict the ground state structure of iron, although the error 
is quantitatively small @]. 

In an attempt to improve upon LDA, a dependence of the exchange-correlation energy 
on the derivatives of the electronic density can be introduced. A simple and systematic 
extension of LDA is the gradient expansion approximation (GEA) H], 



GEA 



[p] = I p(r)6L°^(p(r))dr + J 5.e(p(r))|Vp(r)|Mr + • • • , (2) 

which is is asymptotically valid for densities that vary slowly over space. For finite systems, 
the terms in the expansion of the exchange energy of order greater or equal than four 
diverge while the exchange potential diverges already at second order. By using a constant 
Bxc determined variationally for each atom, Herman et al. showed that the gradient 
expansion up to second order yields improved total energies of atoms |p. An improvement 
in the exchange energies is also obtained when the exact expansion coefficient for exchange, 
-Bx, is used in an exchange-only calculation |^. However, the inclusion of correlation in 
the second order expansion gives energies that are less accurate than LDA 

The analysis of this failure has lead some authors PJTO| to introduce a generalization 
of the gradient expansion, whose generic functional form (here restricted to second-order 
derivative) is 

^x^c^^ [P] = / P(r) et^\p{r), |Vp(r)| , VV(r)) dr. (3) 
Many different GGA's have been proposed since then P,p]T|-[l8| but none of them is clearly 



better than the others. Other classes of approximate functionals have also been tested. 



including schemes that remove the self- interaction of electrons [19,20|, the weighted den- 



sity approximation [^,0 and the average density approximation [^] that attempt to 
model the exchange-correlation hole. However, generalized gradient approximations have 
recently received particular attention both in the physics and chemistry communities since 
they do not add appreciable computational complexity to the LDA scheme. 

The GGA functionals seem to give a good description of several finite systems: they 
significantly improve total energies of atoms [^, as well as binding energies [p^-p6[] and 
vibrational frequencies |Q of ffist and second-row molecules. They also give better esti- 



mates of bond-lengths and binding energies of weakly bound systems such as HA and JIB 



homonuclear dimers ||2^ and hydrogen bonded systems [Q, which are greatly overbound 
in LDA. Lattice constants, bulk moduli and cohesive energies of simple metals and 
transition metals [E^ are also improved. For iron, the correct ferromagnetic bcc ground 



state is predicted [^. On the other hand, the lattice constants and the bulk moduli of 
the semiconductors Si, Ge and GaAs are less accurately predicted by the GGA's than 
LDA, the effect being particularly large for the bulk moduli Calculations for Ad 



and 5c? metals indicate that, while LDA tends to slightly underestimate lattice constants 
and overestimate bulk moduli, GGA's often overcorrect, sometimes leading to results in 
worse agreement with experiment |^|. The tendency of the GGA's to yield smaller bind- 



ing energies and larger bond-lengths than LDA is explained by the fact that the GGA's 
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favor inhomogeneity in the density. Since isolated atoms are more inhomogeneous than 
molecules or solids, the energy of the atoms is lowered more, resulting in smaller bind- 
ing energies. Similarly, for weakly bound systems, larger bond-lengths lead to increased 
inhomogeneity and are therefore favored by the GGA's. 

Studies that compare the predicted structural and energetic properties of various sys- 
tems with the corresponding experimental quantities add very little to the understand- 
ing of the reasons for the successes or the failures of the approximate functionals. In 
order to gain more insight into the limitations of the approximate functionals, many 
researchers |]32|-^ have compared accurate density functional properties, calculated by 
more sophisticated and more computationally demanding methods, with the correspond- 
ing density functional properties calculated from the approximate functionals. In this 
review, we limit ourselves to some of our earlier work and discuss four systems, for which 
it is possible to obtain exact or very accurate wave functions, as illustrative examples: 
a model system of two electrons in a harmonic potential and the He, Be and Ne atoms. 
While, we mostly concentrate on the description of our previous work, there is a new 
ingredient in this paper: the separation of the exchange-correlation potential of many- 
electron systems into separate exchange and correlation components according to the DFT 
definition. 

The outline of the rest of the paper is as follows. In Sec. 0, we briefly introduce density 
functional theory and its Kohn-Sham formulation. In Sec. ^ we list some properties of 
the exact exchange-correlation functional and we determine which properties are satisfied 
or violated by LDA and various GGA's. The short and long distance asymptotic behavior 
of the exchange-correlation potential and local exchange energy per electron is also dis- 
cussed. In Sec. we derive the formulae used to determine the exact density functional 
quantities and, in Sec. H, we describe the systems we study and the corresponding densi- 
ties. The exact exchange and correlation potentials, local exchange energies per electron 
and components of the total energy of these systems are compared with the corresponding 
approximate LDA and GGA quantities in Sees. ||, ^ and ||, respectively. In Sec. ^ we 
summarize our conclusions and discuss some prospects for better approximate functionals. 
In Appendix ^ we list the approximate functionals used in the comparison. 



2. THEORETICAL BACKGROUND 

Density functional theory provides an expression for the ground state energy of a sys- 
tem of interacting electrons in an external potential as a functional of the ground state 
electronic density . Let us assume for simplicity that the spin polarization of the system 
of interest is identically zero. In the Kohn-Sham formulation of density functional theory 
P], the ground state density is written in terms of single-particle orbitals obeying the 
equations in atomic units {fi = e = m = 1): 



2V' + t;ext(r) + 



/ 




dr' + t;xc ([p] ;r) \ = eiip, 



(4) 
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where 

N 



P(r) = El^^Wr- (5) 

The electronic density is constructed by summing over the lowest energy orbitals where 
A'^is the number of electrons, fext(r) is the external potential. The exchange-correlation 
potential v^c {[p] ] r) is the functional derivative of the exchange-correlation energy E^c [p] 
that enters in the expression for the total energy of the system: 

E = -\flJ V^^. dr + I p (r) v,^, (r) dr + 1 II y^dr dr' + E^, [p] . (6) 

The exchange-correlation functional is written as the sum of two separate contributions 
for exchange and correlation, 

^xc [p] = [p] + [p] . (7) 

The definition of the exchange energy is in terms of the non- interacting wave function $0) 
the Slater determinant constructed from the Kohn-Sham orbitals, as 

[p] = ($o| ^ee I <^o)-Uf Ppp^dr dr', (8) 



r — r 



where V^e is the electron-electron interaction. This definition differs from the conventional 
quantum chemistry definition of E^ as the exchange energy in a Hartree-Fock calculation, 
given by the same expression as in Eq. ^ but with the Kohn-Sham determinant replaced 
by the Hartree-Fock determinant. The separation of the exchange-correlation energy 
functional into the separate exchange and correlation components yields a corresponding 
splitting of the exchange-correlation potential into {[p] ; r) and Vc {[p] ; r). In this for- 
mulation, the essential unknown quantity is the exchange-correlation energy E^c [p]- If 
the functional form of i?xc [p], and consequently the exchange-correlation potential, were 
available, we could solve the A''-electron problem by finding the solution of a set of single- 
particle equations. However, as mentioned in Sec. |I], the exact functional form of E^^ [p] 
is not known and it is necessary to make approximations for this term. 

3. PROPERTIES OF THE EXACT DENSITY FUNCTIONAL 

In order to understand if it is possible to construct an approximate functional that 
does not suffer from the weaknesses of the existing GGA's, it is useful to study the known 
properties of the exact density functional and which of these are violated by the commonly 
used approximate functionals. There are actually a surprisingly large number of known 
properties p8|-[53|] and some of these are listed in Table |1|. The GGA functionals compared 



in the table can be found in Appendix ^ 

The first group of properties in Table |I| are sum rules that are satisfied by the exchange 
and correlation holes. The exchange-correlation energy functionals can be interpreted as 
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Table 1 

Some properties of the exact exchange- correlation functional. Comparison with approxi- 



mate functionals. Taken from Ref. and modified 





Property 


E 


LDA pLM 
xc -'-'xc 


pPW91 

xc 


pB88 

X 




LYP 

c 








1 

1 


15 


1 H 




11] 


[18 


fin 1 


171 


1 


/Ox(r,r') <0 




Y 




Y 










2 




,r')dr' = -1 




Y 




Y 










3 


/Pc(r 


,r')dr' = 




Y 




Y 










4 


[p] < 




Y 


Y 


Y 


Y 


Y 






5 


Ec[p] 


< 




Y 


N 


N 






N 


N 


6 


eAp] 


,i?xc [p] > -c///3dr 


a 


Y 


N 


Y 


N 


N 






7 


E^ [px] = XE, [p] 




Y 


Y 


Y 


Y 


Y 






8 


E, [px] <XE,[p], X<1 b 




Y 


N 


Y 






N 


N 


9 




ooEc [px] > -OO 




N 


Y 


Y " 






Y 


Y 


10 


liniA- 


jEc [pa] > -oo 




Y 


N 


Y 






Y 


Y 


11 


liniA- 


oo-Ex [pa] > -oo 




N 


N 


Y 


N 


N 






12 




oE^ [pfl > -oo 




V 


IN 


Y 


Y 


V 










oo i^x [Pa^] > -oo 




Y 


N 


y 


Y 


Y 






14 




jE^ [pfxl > -oo 




N 


N 


Y 


N 


N 






15 


liniA- 


oo XEc [pa] > -oo 




N 


yc 


Y 






N 


N 


16 


limA- 


oiE,[pl]=0 




N 


N 


Y 






N 


N 


17 


liniA- 


oo Ec[pZ]=o 




N 


N 


Y 






N 


N 


18 


liniA- 


i^^c [paa] > 




N 


yc 


Y 






yc 


N 


19 




— > — l/2r, r — > oo 




N 


N 


N 


YN 


N 






20 


Ux(r) 


—1/r, r^oo 




N 


N 


N 


N 


N 






21 




Vc{r) finite, r — > 




Y 


N 


N 


N 


N 


N 


N 


22 


LDA limit for constant p 




Y 


N 


Y 


Y 


Y 


N 


N 


23 


GEA limit for slowly varying 


! P 




N 


Y 


N 


N 


N 


N 



" 1.44 < c < 1.68 

^ Note that E^ [pa] < XE^ [p] , A < 1 is equivalent to E^ \p\\ > XEc [p] , A > 1. 
^ But it diverges to +oo. The PW91 GGA can be modified to satisfy this relation [49[. 
"Y" for exponential p(r) but "N" in general, e.g. e^^^{r) —>■ —1/r for a Gaussian density. 
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the energy arising form the interaction of an electron at r and its exchange-correlation 
hole at r', 

B„|,].l||P(£lP^drdr'. (9) 

The exchange-correlation hole is then separated into exchange and correlation, Pxlr, r') 
and Pc{t^, r') respectively, where the exchange contribution comes from the non- interacting 
system according to Eq. ||. The exchange-correlation hole can be expressed as an integral 
over the coupling constant of the density- density correlation function of the interacting 
system |48| . 

The second group of properties provide bounds on the functionals. Condition 6 is 
known as the Lieb-Oxford bound and a tighter version of it is given in Ref. |^ : 

hm l^xe [pa] = inf (vl> |Ke| ^) - ^ // ^^^drdr' > -c / p(r)^/Mr, (10) 

where the density px is obtained by uniformly scaling the density p in all three spatial 
directions: 

PA(r) = AV(Ar). (11) 

This scaled density integrates to the same number of electrons as the unsealed one. 
The scaling of the exchange and correlation functionals when the charge density is 



scaled uniformly ||50| , |5T|| is described by the third group of properties. It now becomes 
evident why it is useful to separate the exchange component from the entire functional 
according to Eq. ^ For exchange, an exact relation (condition 7) exists under uniform 
density scaling that determines how derivatives of the density combine with the density 
in an exchange GGA functional: 

^X^^[P] = / p(r)^/'F(|Vp(r)|Mr)^/^VV(r)Mr)^/^...)dr. (12) 

The fourth set of properties consists of the relations under non-uniform scaling of the 
density in one or two of the spatial directions: 

PA(r) = Ap(Ax,?/,2;), pZ,{r) = X"^ p {Xx, Xy, z) . (13) 

The non-uniform scaling relations for the correlation energy are here given in the tightest 
form and are derived under the assumption of the existence of a Taylor series in A or 
1/A [^]. A summary of the scaling relations and their derivation is provided in Ref. . 

The fifth group of properties describes the long-distance asymptotic behavior (for finite 
systems) of the exchange potential and the local exchange energy per electron (see Sec. ||) 
and the short- distance behavior of the exchange and correlation potential. 

The last two conditions are that the functionals must reduce to the LDA functional in 
the limit of a homogeneous density and to the correct second-order expansion for a slowly 
varying density. 
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We also mention an additional relation not included in the table, the convexity con- 
straint which is given by 



d'^B [p + eAp] 



+ 



e=0 



Ap(r)Ap(r- 



drdr' > 0, 



r — r 



(14) 



for arbitrary Ap integrating to zero, where 



S[p] = lim-E.e[pA] 



(15) 



This constraint is very stringent: it is violated by the LDA, by the Perdew-Wang '91 
GGA and, likely, by any other GGA functional |^ . 

From Table |^, it appears that most of the functionals violate the asymptotic behavior 
of the exchange-correlation potential and local exchange energy per electron. Within 
GGA, the exchange-correlation potential is given by the functional derivative of E^^^ [p] 
(Eq. I): 



dp 



- V- 



dV'p^ 



(16) 



p{r),Vp(r),V2p(r) 



where e^c is the exchange-correlation energy density, e^c = P^xc- As pointed out in our 
earlier papers [Q,^], any GGA that has terms containing Vp but no higher derivatives of 
p must yield an exchange-correlation potential that diverges at nuclei. If p ~ constant + 



near the origin, then |Vp| 



„s-l 



. If the leading behavior of Cxc in the gradient, near an 



extremum of the density, is iVpp, then 
9e« 



dVp 
and 

V 



iVpf 



~ ^is-l)im-l) 



(17) 



(s-l){m-l)-l 



Therefore, the exchange-correlation potential diverges at the origin if 
(s - l)(m- 1) - 1 < 0, 



(19) 



where s and m are positive. Thus, the exchange-correlation potential diverges at extrema 
of the density for all values of s if m < 1, as is the case for the Wilson-Levy potential []TB||, 



and for all values of m if s < 1. At nuclei, s = 1 so that the exchange-correlation potential 
always diverges. 

As far as the large distance asymptotics are concerned, any GGA that includes no 
higher than first derivatives of the density cannot simultaneously satisfy both the correct 
— 1/r behavior of the exchange-correlation potential and the correct — l/2r behavior of 
the local exchange- correlation energy per electron |1^ . 
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It is however possible, as mentioned in Refs. ||45|j55 
density in an appropriate way 



to include the Laplacian of the 
and construct a GGA that has the triple advantage that 
it satisfies both of the long-distance asymptotic conditions and also does not suffer from 
a spurious divergence at the nuclei. There are an infinite number of possible variations 
in the functional form that preserve the desired short and long range asymptotic behav- 
iors. We are presently searching for the most physically reasonable functional expression 
and for a universal set of parameter values in such a functional. In the same line of re- 
search, Jemmer and Knowles have recently published a simple expression for an exchange 
functional depending on the Laplacian of the density and reproducing the correct long 
and short distance asymptotics of the exchange-correlation potential and local exchange- 
correlation energy per electron However, they conclude that such a functional is not 



a suitable choice as a general purpose density functional. The importance of Laplacian 
terms has also been pointed out by Engel and Vosko [^] who show that inclusion of 
Laplacian terms in a fourth order Taylor expansion of the exchange- correlation functional 
results in improved exchange potentials for atoms and jellium spheres. However, their 
expression does not obey any of the above three asymptotic conditions. 



4. DETERMINATION OF ACCURATE DFT QUANTITIES 

For atomic systems, it is possible to determine accurate exchange-only quantities by 
using the solution of the optimized effective potential (OEP) method which represents 
the exact solution in an exchange-only DFT scheme |2^. However, these quantities are 
determined not for the exact density but for the OEP density which corresponds to the 
self-consistent solution in an exchange-only approach. The performance of exchange- 
only approximate functionals and the corresponding exchange potentials can therefore be 
routinely checked for atomic systems by using the OEP method. Approximate quantities 
in exchange-only GGA, evaluated either for the OEP density or the self-consistent density, 
can be compared with the corresponding OEP quantities. 

Obtaining reliable exchange-correlation potentials and energies is instead a more dif- 
ficult task. One has to generate an accurate density and then compute an exchange- 
correlation potential that yields the desired density as a solution of the Kohn-Sham equa- 
tions (Eq. 1^ and In this context, researchers have used charge densities, of varying 



degrees of accuracy, calculated by quantum chemistry methods for atoms |32-41| and 



molecules as well as Quantum Monte Carlo methods for atoms P5| , p7| and for a 



model semiconductor [^. The subsequent inverse problem, namely the search of the cor- 
responding exchange-correlation potential, has been performed using a variety of different 
techniques. For example, in the special case of the singlet ground state of a two-electron 
system, the exchange-correlation potential can be obtained simply from Eq. ^ while, for 
systems with more than two electrons, v^c can be determined by expanding it in a com- 
plete set of basis functions and varying the expansion coefficients such that Eqs. |^ and || 
yield the accurate density[0,^. Whatever the technique used, an accurate density is a 
key ingredient for the determination of the exchange-correlation potential, because small 
errors in the density are greatly magnified by the inversion procedures used to obtain the 
potential. 
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By knowing the exchange-correlation potential, we are able to calculate the Kohn-Sham 
orbitals and, if the total energy of the system can be estimated, the exchange-correlation 
energy is obtained by inversion of the expression for the total energy (Eq. 

E.c[p]=E-^J2j V^^VV. ^^^^dr dr' - / p (r) t;ext(r)dr. (20) 

i=l r- ^ \ 

Since the Kohn-sham orbitals are known, the exchange energy (Eq. ^ can be calculated 
as 

[p] = -7:1^1^ Sm,,,m.^ / / ^r— — ^ dr dr , (21) 

^ i=i j=i |r - r I 

where the 5-function is over the spin quantum numbers of the i-th and j-th spin-orbitals. 
The correlation energy is finally obtained as the difference of the exchange-correlation 
energy (Eq. ^0|) and the exchange energy (Eq. |2T|) : 

[p] = ^xc [p] - E, [p] . (22) 

The functional derivative of the exchange energy (Eq. ^) with respect to the density 
yields the exchange potential. By knowing the exchange-correlation and the exchange 
potentials, we obtain the correlation potential as the difference 

{[p] ; r) = t;xc ([p] ; r) - v^{[p] ; r) . (23) 



We show in Sub-sec. ^]T] that it is easy to calculate the exchange potential for a system 
of two electron in a singlet state. The procedure is instead much more involved for a 
many-electron system, as we outline in Sub-sec. [4.2| . 

As already introduced in Sec. |l], the exchange-correlation functional can be written as 
an integral over the exchange- correlation energy density, exc(r), or the local exchange- 
correlation energy per electron, e^c (r) : 

E.c [p] = J fixe (r) dr = I p (r) e,, (r) dr. (24) 

The separation of E^^ into exchange and correlation (Eq. |^) corresponds to an equivalent 
decomposition of e^c into and ec- The definition of exc and its components is however 
not unique since two exc's whose difference is a function / with / p (r) / (r) dr = yield 
the same exchange-correlation energy functional. 

We observe that a natural definition of the local exchange energy per electron follows 
from Eq. 

""^"■^ " "2 ^ ^ V'. (r)^,(r) / dr . (25) 



i=i j=i 



This definition results in an that at large distances goes as — l/2r A procedure for 
constructing ec has been recently proposed using first and second order density matrices 
from correlated wave functions and applied to He and H2 ||60[| . 
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In the remaining sections, the accurate exchange-correlation energies and potentials, 
the separate exchange and correlation components and the local exchange energy per 
electron are used to test the accuracy of the approximate exchange-correlation energy 
functionals listed in Appendix |^. 

4.1. Two-electron systems 

For two electrons of opposite spin, there is a single spatial Kohn-Sham orbital ip and is 
simply related to the electronic density (Eq. |^) as 



p(r) 



1/2 



(26) 



Knowing the exact density, the exchange- correlation potential can be obtained by inver- 
sion of the Kohn-Sham equation (Eq. 

vU[p] ; r) = exs + - v^r) - / A^dr'. (27) 

2 ip J |r — r I 

The eigenvalue eKs is equal to minus the ionization energy if we impose that the exchange- 
correlation potential goes to zero at infinity. Therefore, the exchange-correlation potential 
is completely determined. From Eq. it follows that the exchange energy is given by 

and the functional derivative of this expression with respect to the density yields the exact 
exchange potential: 

-x([p];r) = -i/|^dr'. (29) 

Observe that fx {[p] ', r) simply follows from the condition that it cancels the self- interaction 
term in the Hartree potential. 

For a system of two interacting electrons in a singlet state, Eq. for the local exchange 
energy per electron reduces to 



1 f p(r' 



-dr. (30) 



The same expression could have also been obtained from Eq. |2^. 

4.2. Many-electron systems 

While the separation of the exchange-correlation potential into exchange plus cor- 
relation is quite simple in the two-electron systems, this decomposition according to 
the DFT definition has never been obtained for many-electron systems. In previous 
work [^^^^^^ , the exchange potential was defined as the difference of the effective 



Kohn-Sham potential yielding the Hartree-Fock density and the sum of the Hartree and 
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the external potentials. The correlation potential was then obtained as the difference 
of the exchange-correlation potential corresponding to the exact density and the above 
potential. Note that this separation does not correspond to the DFT definition of Eq. |23|: 
it involves two densities, the exact and the Hartre-Fock densities, while the potential used 
for exchange is only approximately equal to the exchange potential corresponding to the 
Hartree-Fock density (although very close to it) since it is not the functional derivative 
with respect to the density of the exchange energy evaluated for the orbitals obtained 
from the effective potential yielding the Hartree-Fock density. Here, we obtain the correct 
separation of the exchange and correlation components, according to their DFT defini- 
tion, for the exact density. The difference between the exchange potential we obtain and 
the approximate "exchange" potential described above is almost not detectable on the 
scale of the exchange potential. On the other hand, on the more expanded scale of the 
correlation potential, the difference is clearly visible, although the shapes are very similar. 
We follow Gorling and Levy in showing how to separate the exchange-correlation 



potential into exchange and correlation. We consider a spin unpolarized system. If we 
assume that the density p is non-interacting 7>representable, it can be expressed as in 
Eq. ^ in terms of single-particle orbitals {ipi} of the Kohn-Sham potential (r), 

V, (r) = ^ext (r) + j y^f^^' + ^xc ([p] ; r) . (31) 

We evaluate the functional derivative of the exchange energy functional with respect to 
the Kohn-Sham potential as 

On the other hand, since the exchange functional can be written as a function of the 
orbitals (Eq. ^1]), we also have 

'^^x [p] ^ f [ (5E^5U^ 5E^[p\ 5ij:{v') \ , 



If we combine Eqs. p2 and 33 , we obtain the integral equation 



t;.([p];r')/C(r',r)dr'= Q(r), (34) 

where the kernel /C (r', r) and the right hand side Q (r) depend on the orbital {ipi} and 
their functional derivative with respect to the potential Vg (r). This integral equation is 
equivalent to the one solved in the OEP method with the KS orbitals replaced by the OEP 



orbitals |^0[- The functional derivatives of the orbitals Sipiiy) / Sv^ir') can be expressed in 



terms of the Green's function Gi (r, r') as 

^ = -G,(r,r')^.(r'), (35) 
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where Gi (r, r') satisfies the differential equation 

i V2 + ^;,(r) - e,) (r, r') = 5 (r - r') - U^Wd^')- (36) 



By knowing the exchange-correlation potential, the KS orbitals and eigenvalues, we can 
compute the Green's functions {Gi} and, consequently, the kernel /C and the function Q. 
If we express the exchange potential as a linear combination of basis functions, Eq. ^ 
can be rewritten as a non-homogeneous set of linear equations for the coefficients of the 
expansion of the potential in the basis set. More details on the separation procedure can 
be found in Ref. [^. 

5. SYSTEMS STUDIED AND CORRESPONDING DENSITIES 

Experimental and computed energies are often compared in order to assess the accu- 
racy of approximate density functionals. However, this procedure is not very reliable 
since energies are integrated quantities and are therefore subject to cancellation of er- 
rors. In order to better understand the performance of approximate exchange-correlation 
functionals, we will examine not only total energies but also the components of the en- 
ergy, the self-consistent charge densities, the exchange-correlation potentials and the local 
exchange energies per electron. We will compare these quantities, obtained from approx- 
imate schemes, with the corresponding accurate quantities derived from (1) the exact 
solution of a model system of two interacting electrons in a harmonic potential (we call 
this system "Harmonium"), (2) a nearly exact wave function for the He atom and (3) 
accurate quantum Monte Carlo calculations for the Be and Ne atoms. 

5.1. Harmonium 

A simple two-electron system where the exchange-correlation potential and energy is 
calculable exactly consists of two interacting electrons in a harmonic potential. The model 
can be solved analytically for a discrete and infinite set of values of the spring constant 
k [^] and the degree of correlation within the system can be varied by simply tuning the 



value of k. 

In Fig. we show the density in a low correlation case, corresponding to a spring 
constant k = 0.25 a.u., and in a high correlation case with a spring constant k ^ 3.6 x 10^^ 
a.u. In the low correlation case, the density has a maximum at the origin while, in the high 
correlation case, the density has a local minimum at the origin and an absolute maximum 
at a finite distance from the origin Although the high correlation density is too 



low to be physically relevant for electronic structure calculations ((r^) = 57.9 a.u.), the 
system is interesting because it is strongly correlated and its density differs qualitatively 
from an atomic density due to the presence of a maximum at a finite distance from the 
origin. Consequently, in the present paper, we will only discuss the high correlation limit. 
For an analysis of the system with k = 0.25 a.u., see Ref. p^. 



5.2. He atom 

It is necessary to employ very accurate wave functions in order to observe the true 
short-distance and long-distance behavior to 0{l/r) of the exchange-correlation potential. 
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Figure 1. Ground state electronic densities: the spring constants are k = 0.25 a.u. (left, 



low correlation case) and k 
Ref. 1 5^] and modified. 



3.6 X 10 a.u. (right, high correlation case). Taken from 
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Figure 2. Comparison of the LDA, GGA (PW91), HF and exact densities of He (left). 
Errors in the LDA, GGA, and HF densities of He (right). Taken from Ref. and 
modified. 
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For the He iso-electronic series, it is possible to achieve this accuracy. The form of the 
wave function is described in Refs. The energy for He corresponding to this 

wave function was estimated to be about 1 part in 10^^. With this nearly exact wave 
function, the exchange-correlation potential is resolved with unprecedenced accuracy: it 
is possible to observe not only the — 1/r behavior of the exchange potential but also 
the asymptotic behavior of the correlation potential, — 9/(4Z^r'^), predicted by Almbladh 
and von Barth |65 |. 

In Fig. 1^, we compare the density obtained with our accurate wave function with the 
self-consistent densities obtained from LDA, the Perdew-Wang '91 GGA[Q and Hartree- 
Fock (HF). The LDA density is less peaked than the true density while the HF density is 
very close but slightly more peaked than the true density. We also show the error in the 
LDA, GGA, and HF densities. The GGA density is more accurate than the LDA density 
but less accurate than the HF one. Note that the error in the HF density is considerably 
smaller than the LDA and GGA errors; however, as shown below, this is true only in the 
core region for heavier systems. 



5.3. Be atom 

For many-electron systems, an accurate determination of the charge densities is ob- 
tained by combining the results from variational Monte Carlo (VMC) and diffusion Monte 
Carlo (DMC) methods In VMC, the square of a trial wave function is sampled and 



the expectation values for the trial wave function are calculated as averages over the sam- 
pled configurations. An advantage of this technique versus other conventional quantum 
chemistry methods is that the many- dimensional integrals can be performed for any given 
form of the trial wave function. Therefore, one is freed from the constraint of having to 
express the trial wave function in some restricted form such as a linear combination of 
determinants of single-particle orbitals. The starting point of the DMC method is also a 
good trial wave function from which DMC projects out an improved estimate of the true 
wave function. 

To calculate the charge density, the integral of the square of the true wave function 
over all but one of the electrons has to be estimated. VMC samples the square of the trial 
wave function. DMC samples the product of the trial wave function and the true wave 
function. If the error in the wave function is of 0{e), the errors in the densities obtained 
from VMC and DMC are of the same order. A density with an error of O(e^) (that we 
will refer to as the quantum Monte Carlo (QMC) density) is obtained by taking twice the 
DMC density and subtracting the VMC density ||67[|. The QMC density is obtained in the 
form of an histogram and then fitted to a sum of products of monomials, exponentials 
and an appropriate asymptotic function [|6^] . 

In Fig. ^, we show the densities (left) obtained from LDA, the Perdew-Wang '91 GGA, 
HF and QMC for Be and the errors in these densities (right), with the QMC density as 
the reference. The first thing to notice is that all densities are very similar. The LDA 
and HF errors are of comparable magnitude in the valence region whereas, as in the 
case of He (Fig. the HF error is smaller than the LDA error in the region near the 
nucleus. A hand- waving explanation of this behavior is that, in the core region, exchange 
dominates correlation and HF, by definition, is exact for exchange only. The Perdew- 
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Figure 3. Comparison of the LDA, GGA (PW91), HF and "exact" (QMC) densities of 



Be (left). Errors in the LDA, GGA, and HF densities of Be (right). Taken from Ref. [^ |. 



Wang '91 density is somewhat more accurate than the LDA density. The GGA potential 
has a spurious negative divergence at the nucleus that is partially responsible for the 
improvement in the density. The LDA density is too low at the nucleus and the negative 
divergence of the GGA potential results in increased charge density there, a feature we 
have observed for Ne as well. 



6. A CLOSER LOOK AT EXISTING APPROXIMATIONS: V^c, AND Vc 

The accurate exchange-correlation potentials and their separation in exchange and cor- 
relation components were obtained for the four systems described in this paper, using 
the techniques mentioned in Sec. The approximate exchange- correlation potentials are 
evaluated numerically as functional derivatives of the approximate exchange-correlation 
energy functional (Eq. |l^) and compared with the corresponding accurate quantities. 

6.1. Harmonium 

In Fig. H, we plot the exact exchange potential and several approximate exchange po- 
tentials evaluated for the exact density of the model system. The approximate potentials 
differ from the exact potential over all the spatial range. At large r, all the potentials 
do not reproduce the correct — 1/r asymptotic behavior: they go exponentially to zero 
with the exception of the Becke functional which behaves as (— constant/r^) |T^ and of 
the Langreth-Mehl functional which diverges. At the origin, the exact exchange potential 
has a minimum because it is simply proportional to the Hartree potential through Eq. 
while the approximate exchange potentials exhibit a maximum (they do not diverge since, 
from Eq. [T^, s = 2 and m > 1). The Perdew-Wang '91 exchange potential is very close 
to the Becke '88 potential except that it displays additional oscillations at the extrema 
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Figure 4. Exact and approximate exchange (left) and correlation (right) potentials for 
the spring constant k ~ 3.6 x 10~^ a.u. Taken from Ref. 0| and modified. 



of the densities and in the tail of the density. In addition, as mentioned earlier, the long- 
distance asymptotic behaviors are different. As pointed out by Engel and Vosko [^], the 
oscillations of the Perdew-Wang '91 potential at the quadratic extrema of the density are 
due to the fact that the effective coefficient of the |Vpp term increases very rapidly from 
the known value at ^ = 0, specified by the second-order gradient expansion, to a 2.1 times 
larger value at only ^ = 0.04, where ^ = [Vp/(2fcp'p)]^ and kp is the Fermi wave vector. 
Observe that the Becke '88 potential does not have spurious oscillations at the extremum 
of the density nor in the tail of density. 

In Fig. ^, we also show the exact correlation potential and several approximate corre- 
lation potentials evaluated for the exact density. At large r, the approximate potentials 
go to zero exponentially with the exception of the Wilson-Levy potential which goes to 
a positive constant and the Langreth-Mehl functional which diverges. As predicted by 
Eq. [l^, the Wilson-Levy potential diverges at the origin: the energy density, exc(p, |Vp|), 
depends on |Vp| as (9(|Vp|)) for small values of the gradient of the density (Eq. |19|, 
m = 1). All the other functionals assume a finite value at the origin (Eq. s = 2 and 
m > 1). The Wilson-Levy functional has a discontinuity at the maximum of the electronic 
density: the functional derivative of the exchange-correlation energy contains the sign of 
the radial derivative of the density and is discontinuous if the energy density depends on 
the absolute value of gradient of the density as (9(|Vp|)), for small values of |Vp|. As 
in the case of exchange the Perdew-Wang '91 correlation potential shows an oscillatory 
behavior near the maximum at a finite distance from the origin and also in the tail of the 
potential. 
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6.2. He atom 

For atomic systems, the exchange potential is the dominant part of the exchange- 
correlation potential. It shares with it most of its characteristics: the asymptotic behavior, 
the intershell bump (see Sub-sec. p. 3D and the finite short-range behavior. 

In Fig. I we compare the exchange potential obtained from the LDA and the various 
GGA's evaluated for the exact charge density with the exact potential obtained from 
Eq. The approximate potentials differ significantly from the exact one. At large 
r, all the potentials do not reproduce the correct — 1/r asymptotic behavior: they go 
exponentially to zero with the exception of the Becke '88 functional which behaves as 
(— constant/r^) [0 and of the Langreth-Mehl functional which diverges. The Perdew- 
Wang '91 functional has a spurious minimum between 2.5 and 3.5 oq. 

At the origin, the exact exchange potential has a quadratic minimum because it is 
simply proportional to minus the Hartree potential. In contrast all the proposed GGA 
exchange potentials diverge at the nucleus while the LDA exchange potential has a finite 
value and slope there. Hence, very close to the nucleus, the various GGA potentials are an 
even poorer approximation to the true potential than LDA but the more negative values 
of the GGA potentials at short and intermediate distances is a step in the right direction. 
As shown in Eq. ^ any approximate functional that has terms containing the gradient 
but no higher derivatives of p must yield an exchange-correlation potential that diverges 
at nuclei since, at nuclei, s = 1. 

In Fig. 1^, we also compare approximate correlation potentials of He to the exact one. In 
the case of an n-electron atom or ion with orbitally non-degenerate n and n — 1 electron 
ground states, the correlation potential goes as — a/2r'^ at large distances, where a is the 
dipole-polarizability of the n — 1 electron system p5[]. In the case of the He iso-electronic 
series, the dipole-polarizability of the residual 1-electron system is a = 9/(2Z'^). Due 
to the rapid fall-off of 1/r^, this asymptotic behavior is not discernible in Fig. |] but is 



evident in a plot of r^Vc |55|. At large distances, the approximate potentials go to zero 
exponentially with the exception of the Wilson-Levy potential which goes to a positive 
constant and the Langreth-Mehl potential which diverges. All the GGA correlation po- 
tentials diverge at the origin as follows from Eq. |19|. Of the GGA's proposed so far, 
only the Lee- Yang-Parr GGA contains Laplacian terms but not in the form necessary to 
eliminate the divergence at the nucleus. 



6.3. Be and Ne atoms 

In Fig. we show the accurate exchange potentials of Be and Ne and compare them 
with the approximate exchange potentials from the local density approximation and sev- 
eral generalized gradient approximations, all of them being evaluated for the correspond- 
ing accurate QMC density. 

At short and large distances, all GGA's are similar and qualitatively wrong and the same 
analysis applies as in the case of the He atom. The Perdew-Wang '91 exchange potential 
has a spurious dip at about 6 a.u. (outside the plotted range) for Be and at about 3.3 a.u. 
for Ne. This dip is not present in either Perdew-Wang '86 or the Becke '88 potentials. 
The exact exchange potential has a peak in the intershell region (around 1 a.u. for Be 
and around 0.3 a.u. for Ne). There, the various GGA potentials exhibit some differences 
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Figure 5. Exact and approximate exchange (left) and correlation (right) potentials for 
He. Taken from Ref. |^5| and modified. 
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Figure 6. Exchange potentials of Be (left) and Ne (right) from QMC (labeled exact - see 
text), LDA and various GGA's. Taken from Ref. 




Figure 7. Correlation potentials of Be (left) and Ne (right) from QMC (labeled exact - 
see text), LDA and various GGA's. Taken from Ref. 



but they all have a peak at approximately the right position. The GGA potentials are 
more accurate than the LDA potential at distances smaller than the intershell radius, 
except at very short distances where they have a spurious divergence. To summarize, the 
exchange potentials from the GGA's show some tendency for improvement upon LDA in 
the intershell region and at shorter distances but they have a spurious divergence at the 
nucleus and do not improve the long-distance behavior. 

In Fig. 1^, we show the correlation potentials for Be and Ne from QMC, LDA and 
the different GGA's. The LDA correlation potential is everywhere negative, monotoni- 
cally increasing with the radius (since the density is monotonically decreasing) and rather 
smooth. Except for the part of the potential beyond 1.5 a.u. for Be, it bears no resem- 
blance to the exact correlation potential. The GGA correlation potentials exhibit more 
structure but none of these does better than LDA. Nevertheless, as we will see in Section || 
the GGA's yield more accurate correlation energies than does LDA. The true correlation 
potential is everywhere negative for Be but is predominantly positive for Ne. 



7. LOCAL EXCHANGE ENERGY PER ELECTRON: 6^ 

Using the expression in Eq. ^ for the local exchange energy per electron, ex, we present 
accurate ex's for the two-electron systems and the Be atom. Since ex is not uniquely defined 
(see Sec. the differences of approximate ex's from each other and from the "true" one 
are not necessarily indicative of a failure of the approximate functionals. 

In the left panel of Fig. ^ we plot the exact and approximate local exchange energies 
per electron for Harmonium. At large r, the approximate exchange functionals go to zero 
exponentially with the exception of the Langreth-Mehl functional which diverges and of 
the Becke functional which has the wrong asymptotic behavior — 1/r for Harmonium al- 
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Figure 8. Exact and approximate local exchange energy per electron for Harmonium 
[k ^ 3.6 X 10~^ a.u.) (left) and He (right), evaluated for the respective exact densities. 
Taken from Refs. |Q and |^ and modified. 




Local Exchange Energy 
per Electron ot Beryllium 



-1/2r 

Exact 

LDA 

Langreth-Meh 
Perdew-Wang '86 
Perdew-Wang '91 
Becke '88 



0.0 0.5 1.0 



1.5 2.0 2.5 J.O 3.5 4.0 4.5 5.0 5.5 6.0 

r (an) 




Figure 9. Exact and approximate local exchange energy per electron for Be (left) and Ne 
(right), evaluated for the respective exact densities. 
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though it has the correct asymptotic behavior — l/2r for atoms. The Becke functional was 
constructed to have the correct asymptotic behavior for exponentially decaying density 
but it fails for Gaussian densities (see Table |1|). At the extrema of the density, all the 
functionals go to the LDA value. 

In the right panel of Fig. ^ and in Fig. ^ we compare the exact and approximate local 
exchange energies per electron of the He, Be and Ne atoms. At large r, the approximate 
ex go to zero exponentially with the exception of the Langreth-Mehl functional, which 
diverges, and the Becke functional which has the correct asymptotic behavior — l/2r. 
However, as also noted by Engel et al. , the Becke achieves its asymptotic behavior 
at much larger distances than does the true e^- This is because the correction to the 
leading — l/2r term is exponentially small for the true functional but (9(l/r^) for the 
Becke functional. The Becke exchange energy per electron is however closer to the correct 
one over all the range than any other approximate e^. At intermediate distances, the 
various GGA's represent an improvement upon LDA, yielding a deeper than LDA, but, 
at very short distances, they are significantly lower than the true one. 

8. COMPARISON OF THE COMPONENTS OF THE TOTAL ENERGY 

For the model system of two electrons in a harmonic potential, the He iso-electronic 
series and the Be atom, we present a comparison of the accurate values of E^, and 
i?xc with the ones corresponding to approximate functionals evaluated for the accurate 
densities. In order to assess the change that results from employing the self-consistent 
rather than the exact density, we also perform, for Be, self-consistent calculations with 
LDA and the Perdew-Wang '91 GGA and compare the different components of the total 
energy (Eq. ^ with the accurate values. 

For the model system of two interacting electrons in a harmonic potential, we compare 
in Table ^ the values of the exact exchange and correlation energies with those obtained 
from the approximate functionals evaluated for the exact electronic density. The GGA 
functionals, except for the Langreth-Mehl functional, improves considerably upon LDA. 
The major corrections are achieved for the exchange contribution. The failure of the 
Langreth-Mehl functional is a result of the diverging energy per electron at large r. Most 
of the approximate functionals are constructed to give accurate values of E^^ for atoms 
and it is encouraging that they also improve E^c for this rather different model system. 
However, we can see in Tables ^ and ^ that the improvement upon LDA is much larger for 
atomic systems: for instance, the percentage error in the LDA exchange energy evaluated 
for the electronic density of He is 14% and in the LDA correlation energy, —165%, while 
the percentage errors for the other approximate functionals we list in the appendix range 
from —1 to 1% for exchange and from —3 to 19% for correlation. If we consider the sum 
of the exchange energy and the correlation energy, LDA yields a good result due to a 
cancellation of errors and the GGA's do not improve upon LDA. 

In Tables ^and ^, we compare the approximate i?x, Ec and E^c, evaluated for the accu- 
rate densities of the He iso-electronic series and Be respectively, with the corresponding 
accurate quantities. As is well known, LDA yields values of E^ that are too small in 
absolute magnitude and values of Ec that are too large in absolute magnitude, resulting 
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Table 2 

Model system with spring constant k ^ 3.6 x 10~^ a.u. Exchange-correlation energy for 
the approximate functionals listed in Appendix ^. The functionals are evaluated for the 
exact density. The exact results are also reported. The numbers in parentheses are the 
percentage errors. The total energy of the system is 0.0228 a.u. Energies are in Hartree 
atomic units. Taken from Ref. and modified. 



Functional 


Ex 


Ec 


Exc 




Exact KS 


-0.0195 


-0.0064 


-0.0259 




IDA (Perdew-Wang) 


-0.0174(11%) 


-0.0108(-69%) 


-0.0282(- 


-9%) 


Langreth-Mehl 


-0.0217(-11%) 


-0.0119(-86%) 


-0.0336(- 


-30%) 


Perdew-Wang '86 


-0.0198(-2%) 


-0.0096(-50%) 


-0.0294(- 


-14%) 


Perdew-Wang '91 


-0.0195(0.2%) 


-0.0084(-31%) 


-0.0279(- 


-8%) 


Becke '88 


-0.0201(-3%) 








Wilson-Levy 




-0.0145(-127%) 






Lee- Yang-Parr 




-0.0032(50%) 







in a cancellation of errors for Exc- All the GGA's yield improved exchange, correlation 
and exchange-correlation energies compared to the LDA. 

For the He iso-electronic series, the true correlation energies are nearly constant, whereas 
for the Be iso-electronic series there is a linear component in the nuclear charge Z. It is 
difficult for any GGA functional to mimic the correct Z dependence of both iso-electronic 
series. 

In Table ^, we show various quantities, computed by LDA and the Perdew-Wang '91 
GGA, and the corresponding accurate values for Be. LDA and GGA values were obtained 
through a self-consistent procedure and are therefore evaluated for the LDA and the GGA 
density respectively. The change in Ex and Ec due to self-consistency are not negligible. 
The values of Ex and Ec in Table |^ differ from those in Table |^ by a larger amount 
for LDA than for the Perdew-Wang '91 GGA reflecting our earlier observation that the 
GGA yields somewhat more accurate self-consistent densities than does LDA. The last 
two columns show the errors in LDA and GGA. The LDA components of the total energy 
are all smaller in absolute magnitude than the corresponding accurate values. Both the 
total energy and its components have much smaller errors in GGA than in LDA. In 
Fig. ^ we observed that, near the nucleus, the LDA density is too small while GGA gives 
an increased charge density due to the lower value at short distances and the negative 
divergence at the nucleus of the GGA exchange-correlation potential. This increase in 
the GGA density results in an increase in the magnitude of all the components of the 
total energy. Consequently, the GGA yields improvements not only in Exc but also in 
other components of the total energy, despite, or may be one should say because of, 
the spurious divergence of the exchange-correlation at the nucleus. On the other hand, 
quantities related to the single-particle eigenvalues (shown in the last two rows) are not 
improved in GGA relative to LDA. In particular, the highest lying eigenvalue which, for 
the true density functional, should equal the ionization energy, is not much improved. 
This is to be expected since the long-distance behavior of exchange- correlation potential. 
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Table 3: Comparison of exchange and correlation energies of the He iso-electronic series for LDA and the various GGA's with 
the corresponding exact values. The approximate functionals are evaluated for the exact (not the self-consistent) densities. The 
numbers in parentheses are the percentage errors. All numbers are accurate to all digits presented with the exception of the 
"exact" results for Hg"^^^ which probably have an error in the last digit. Taken from Ref. [55] and modified. 



E. 



Ion 


Exact 


LDA 


LM 


PW86 


PW91 


BBS 


WL 


LYP 






[15] 


[9] 


[12,13] 


[14] 


[11] 


[16] 


[17] 


H- 


-0.380898 


-0.337(12%) 


-0.387(-2%) 


-0.400(-5%) 


-0.393(-3%) 


-0.395(-4%) 


— 


— 


He 


-1.024568 


-0.883(14%) 


-1.014(1%) 


-1.033(-1%) 


-1.016(1%) 


-1.025(-0%) 


— 


— 


Be+2 


-2.276628 


-1.957(14%) 


-2.246(1%) 


-2.279(-0%) 


-2.245(1%) 


-2.265(1%) 






Ne+8 


-6.027484 


-5.173(14%) 


-5.940(1%) 


-6.015(0%) 


-5.927(2%) 


-5.982(1%) 






Hg+78 


-49.777931 


-42.699(14%) 


-49.036(1%) 


-49.605(0%) 


-48.889(2%) 


-49.348(1%) 






H- 


-0.041995 


-0.072(-71%) 


-0.049(-16%) 


-0.038(10%) 


-0.032(24%) 




-0.032(25%) 


-0.030(29%) 


He 


-0.042107 


-0.112(-167%) 


-0.050(-19%) 


-0.044(-4%) 


-0.046(-9%) 




-0.042(1%) 


-0.044(-4%) 


Be+2 


-0.044274 


-0.150(-240%) 


-0.025(43%) 


-0.049(-12%) 


-0.053(-21%) 




-0.045(-2%) 


-0.049(-ll%) 


Ne+8 


-0.045694 


-0.202(-342%) 


0.082(279%) 


-0.084(-83%) 


-0.062(-35%) 




-0.047(-3%) 


-0.050(-10%) 


Hg+78 


-0.046538 


-0.326(-600%) 


1.833(4038%) 


-0.299(-543%) 


-0.081(-73%) 




-0.048(-3%) 


-0.051(-9%) 


H- 


-0.422893 


-0.408(3%) 


-0.435(-3%) 


-0.438(-4%) 


-0.425(-0%) 








He 


-1.066675 


-0.996(7%) 


-1.064(0%) 


-1.077(-1%) 


-1.062(0%) 








Be+2 


-2.320902 


-2.107(9%) 


-2.271(2%) 


-2.328(-0%) 


-2.298(1%) 








Ne+8 


-6.073178 


-5.375(11%) 


-5.858(4%) 


-6.099(-0%) 


-5.989(1%) 








Hg+78 


-49.824470 


-43.025(14%) 


-47.203(5%) 


-49.904(-0%) 


-48.970(2%) 
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Table 4 

Comparison of exchange and correlation energies of Be for LDA and the various GGA's 
with the corresponding accurate values. The approximate functionals are evaluated for the 
exact (not the self-consistent) densities. The numbers in parentheses are the percentage 
errors. Energies are in Hartree atomic units. Taken from Ref. ||47[| . 



Functional 



E^ 



Accurate 

LDA (Perdew-Wang) 
Langreth-Mehl 
Perdew-Wang '86 
Perdew-Wang '91 
Becke '88 
Wilson-Levy 
Lee- Yang-Parr 



-2.67398 

-2.32125(13%) 

-2.60576(3%) 

-2.69050(-l%) 

-2.65434(1%) 

-2.66694(0%) 



-0.09615 

-0.22515(-134%) 
-0.10012(-4%) 
-0.09425(2%) 
-0.09496(1%) 

-0.09572(0%) 
-0.09546(1%) 



-2.77013 

-2.54641(8%) 
-2.70588(2%) 
-2.78475 (-1%) 
-2.74929(1%) 



Table 5 

Comparison of LDA, Perdew-Wang '91 GGA and accurate values for various components 
of the Kohn-Sham total energy, the density-weighted integral of v^c and single-particle 
eigenvalues of Be. The approximate functionals are evaluated for the respective self- 
consistent densities. The accurate values probably have errors in the last digit quoted. 
The accurate value of the total energy is from Ref. Energies are in Hartree atomic 
units. Taken from Ref. [4^ . 



Property 


LDA 


PW91 


Accurate 




Apw91 


E = Ts + Een + Eii + E^c 


-14.44647 


-14.64794 


-14.66736 


0.22089 


0.01942 


T 


14.30879 


14.57461 


14.5942 


-0.2854 


-0.0196 


Een 


-33.35610 


-33.65670 


-33.7099 


0.3538 


0.0532 


En 


7.11487 


7.17261 


7.2185 


-0.1036 


-0.0459 


E^c 


-2.51403 


-2.73845 


-2.77013 


0.25610 


0.03168 


E. 


-2.29033 


-2.64428 


-2.67398 


0.38365 


0.02970 


Ec 


-0.22370 


-0.09417 


-0.09615 


-0.12754 


0.00198 


Een + Ell 


-26.24123 


-26.48410 


-26.4914 


0.2502 


0.0073 


Jdv p(r) i;xc(r) 


-3.30614 


-3.50058 


-4.4595 


1.1533 


0.9589 




-3.85609 


-3.91154 


-4.2265 


0.3704 


0.3149 


^2s 


-0.20577 


-0.20719 


-0.3426 


0.1368 


0.1354 
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which strongly influences the highest eigenvalue, is not appreciably improved. All the 
observations made here for the Be atom are also true for the Ne atom (see Refs. [^] and 



9. PROSPECTS 

Although the exchange and correlation energies calculated from the various GGA's 
(with the exception of the Langreth-Mehl GGA for the model system and for the corre- 
lation energies in the He iso-electronic series) are considerably more accurate than those 
from LDA, the situation is not as clear for the exchange and correlation potentials. Since 
and Ec are integrated quantities, it is possible for GGA's to improve upon them with- 
out making comparable improvements to the exchange and correlation potentials which 
contain more detailed information. 

The various GGA exchange functionals have certain features in common. By favoring 
regions with large density gradients, they correctly lead to more negative total energies 
and reduced binding energies relative to LDA. They yield a lower potential than LDA at 
short and intermediate distances from nuclei and approximately reproduce the intershell 
structure - two steps in the right direction - but they also introduce a spurious divergence 
at the nucleus and fail to improve the long-range behavior. Whereas the various GGA 
exchange functionals have some common features, the GGA correlation functionals are 
very different from each other and none of them yields a correlation potential resembling 
the true one. Correlation is more subtle than exchange as manifested in several ways. 1) 
The exchange potential is everywhere negative, whereas the correlation potential can be 
either negative or positive in different regions of space. 2) For He iso-electronic series, 
the correlation potential varies less over space for heavy ions than for light atoms p5 |. 



This is opposite to the behavior in LDA. 3) The correlation potential has qualitatively the 
same behavior for all members of the helium iso-electronic series (negative at short 
distances and long distance, positive at intermediate distances, see Fig. |) but it has the 
opposite qualitative behavior for the closely related model system of two electrons in a 
harmonic potential (Fig. 

As discussed in Sec. |^, any GGA that includes no higher than first derivatives of the 
density cannot simultaneously satisfy both the correct — 1/r behavior of the exchange- 
correlation potential and the correct — l/2r behavior of exchange energy per electron 
at large distances. Moreover, any such GGA must have a spurious divergence in the 
exchange-correlation potential at nuclei. We mentioned that, by including the Laplacian 
of the density in an appropriate way, it is possible to construct a GGA that has the correct 
long-range asymptotics of the exchange potential and energy per electron and the proper 
short-range behavior of the exchange-correlation potential. Nevertheless, it may turn out 
that even GGA functionals containing Laplacian terms are not flexible enough to achieve 
a major improvement and it may be necessary to consider more general functionals. We 
discuss here some possibilities. 

Since a data-base on accurate exchange-correlation potentials has being built in recent 
years and none of the existing GGA functionals passes the test of the comparison with 
the true potentials, more effort is now being devoted to developing schemes that focus on 
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reproducing the main features of the true exchange-correlation potential. As mentioned 
above, the potential controls the quality of the self-consistent density and is the key 
ingredient in evaluating energy gradients ||7I|] , 



dX J ' dX 

whose accurate estimation is for instance necessary in determining potential energy sur- 
faces. 

Instead of starting from an exchange-correlation energy functional and then deriving the 
corresponding potential, some authors start with an expression for the potential which has 



the correct asymptotic behavior p9| , [73| . While it is easy to construct a potential with the 
correct asymptotics, one encounters the problem of how to obtain the exchange-correlation 
energy functional that corresponds to this potential. In fact, the potential may not be 
the functional derivative of an energy expression. If the approximate exchange potential 
satisfies the uniform scaling property {[px] ; r) = Xv^ ([p] ; Ar), a possibility is to use the 



Levy-Perdew relation to assign a corresponding exchange energy |50|: 



E. = Jv^ ([p]; r) (3p(r) + r • Vp(r)) dr. (38) 

The relation is derived by considering a line integral of the potential along the density 
path px (Eq. |Tl]) and is exact only for a potential that is a functional derivative and has 
the correct scaling property: 

E. [p] = I' dX J dr iipx] ; r) = J ([p]; r) (3p(r) + r • Vp(r)) dr. (39) 

The first equality always holds for a potential that is a functional derivative while the 
second one is valid only if the potential has the correct behavior under uniform scaling. 
Van Leeuwen and Baerends ||7T[ discuss why Eq. ^ is a sensible choice for a potential 
that is not a functional derivative but has the proper scaling property and is constructed 
to mimic the true potential at the given density. Unfortunately, such a relation does not 
exist in the case of exchange and correlation since the scaling properties of the correlation 
potential are not homogeneous [^]. Assigning an exchange-correlation energy to a given 
potential becomes more complicated because we need to know the potential not only for 
the density of interest, as in the case of exchange-only, but along some path in density 
space [[fT |. 



Another attempt to generate a DFT scheme yielding a better description of the ex- 
change-correlation potential has as a starting point the OEP method. The OEP method 
provides the exact solution for the problem of exchange-only DFT but, as shown by 
Norman and Koelling |73], it can be generalized to treat exchange and correlation if 



the correlation energy is also expressed as a functional of the Kohn-Sham orbitals. The 
resulting problem is technically difficult to solve but Krieger, Li and lafrate proposed an 
approximation that greatly simplify the approach [^]. Although the resulting exchange- 



correlation potential is no longer equivalent to the functional derivative with respect to 
the density of the original exchange-correlation energy, it closely corresponds to the exact 
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generalized OEP exchange- correlation potential if the expression 5Ey,c [{''Pi}] 
goes to zero at large distances. Since the true exchange energy is known as a functional 
of the Kohn-Sham orbitals (Eq. ^) and can therefore be employed, the approximate 
exchange potential, evaluated according to the procedure of Ref . , has the correct short 
and long-distance asymptotic behavior although the intershell bump is underestimated. 
An improved approximation that yields a more accurate intershell bump is also discussed 
in Ref. ||7^. This approach is of practical interest since it is easier to construct self- 
interaction free, exchange-correlation functionals that depend on the orbitals instead of the 
density. Grabo and Gross |]75[ recently implemented this scheme with the approximation 
of Krieger, Li and lafrate and the orbital-dependent correlation functional of Colle and 



Salvetti [76 



A different approximate scheme that cures the long-distance asymptotics of the ex- 
change-correlation potential is obtained by combining a given approximate exchange- 



correlation functional with a self-interaction correction (SIC) procedure |T9[. A GGA de- 
pending on the Laplacian of the density could be easily constructed so that the exchange- 
correlation potential does not have a spurious divergence at nuclei and could then be 
implemented in a SIC scheme to yield a potential with also the correct long-range asymp- 
totic behavior. 
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A. APPROXIMATE FUNCTIONAL FORMS OF E^^ [p] 

The second order generalized gradient approximation exchange-correlation energy is 
written as 

^x^c""^ [P] = / dre.e(p(r), | Vp(r)| , VV(r)). (40) 
We define the following density dependent variables: 

We present the spin unpolarized version of the approximate functionals, ( = (pt — /p = 
0. All the parameters that appear in the following functionals are in atomic units. 
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LDA exchange functional 



^LDA 



where = -(3/4) (S/vr)^/^ 

LDA correlation functional (Perdew-Wang [1l5|] ) 



(41) 



.LDA 



-2ap (1 + aiTs) log 



1 + 



1 



3/2 



(42) 



where a = 0.0310907, = 0.21370, A = 7.5957, /^s = 3.5876, /?3 = 1.6382 and 
0.49294. 



RPA correlation functional 1177 



Cc = -pc^/(r,/rP), 
m = (l + .^)ln(l + l) + ^-.^-l. 



where = 0.0252 and = 30. 

Langreth-Mehl exchange-correlation functional [^] 



(43) 



e-A _ J_YpI (I + 18 f 



..RPA 



p4/3 V9 



(44) 
(45) 



where F = b | Vp| /p^/^ b = {97rf^^ f,a = 7r/(16(37r2)4/3) and / = 0.15. 



Perdew-Wang '86 exchange functional [ |12[| 



= e^°^(p) 1 + 0.0864— + 6/ + cs^ 
V m 



(46) 



where m = 1/15, 6 = 14 and c = 0.2. 
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Perdew-Wang '86 correlation functional ||13|| 



4/3 



where 



$ = 1.745/ 



P' 



~Cc(p = oo) |Vp| 



(47) 



Cc(p) P^/6 



C.(p) = C^ + C,r.^C,rl 



1 + Csr, + C^rl + Crr^ 



and / = 0.11, Ci = 0.001667, C2 = 0.002568, C3 
8.723, Ce = 0.472, C7 = 7.389 ■ lO-^. 



0.023266, Ci = 7.389 ■ 10-^ C5 



Perdew-Wang '91 exchange functional ||14|| 



ex = e.^°^(p) 



l + flissinh ^ {a2s) + (^a^ + a^e -^^^^^ 
1 + ai s sinh~^(a2s) + a^s"^ 



where ai = 0.19645, 03 = 7.7956, 03 = 0.2743, = -0.1508 and = 0.004. 
Perdew-Wang '91 correlation functional ||14| 



(4J 



where 
H = 
A = 



2a 
2a 

J 



log 



1 + 



2a t^ + At^ 
1 + At^ + AH^ 



+ C,o [Cc(p)-Cei] fe 



2^-100 



(49) 



and a = 0.09, /3 = 0.0667263212, Cco = 15.7559, C^i = 0.003521. The function C,(p) is 
the same as for the Perdew-Wang '86 correlation functional, e^^'^(p) = e^^'^{p)/p. 



Becke '88 exchange functional ||11|| 



ex = er"(p) 



/3 



X 



2V3^xl + 6/?xsinh-^(x) 
where x = 2 {Qir^f^^ s = 2^3 |Vp| /p^^^ A^ = -(3/4) and (3 = 0.0042. 



(50) 
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Engel-Chevary-Macdonald-Vosko exchange functional [| 



_ LDA. ^l±^li + M! 

where ^ = and ai = 27.8428, 02 = 11.7683, h = 27.5026 and 62 = 5.772^ 
Wilson-Levy correlation functional [Iir 



(51) 



ap + 6 |Vp|/p 



1/3 



' C + f/|Vp|/(p/2)4/3+r,' 

where a = -0.74860, b = 0.06001, c = 3.60073 and d = 0.90000. 



Closed shell Lee- Yang-Parr correlation functional []T7| 



(52) 



en = -a 



where 



l + dp- 



tw — - 



p 



-cp-1/3 



and Cf = 3/10 {Sn'^f^'^, a = 0.04918, b = 0.132, c = 0.2533 and d = 0.349. 
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